Collaborators: Tomo Tanaka, John Eykelenboom.

shiny app

1 Data

In our data we distinguish four states, marked by colour: blue, brown, pink, red. For the purpose of the model we merge blue and brown together into one state. Here is the result of the experiment from 38 cells.

The next figure shows the proportions of each colour as a function of time. The curves are smoothed with a running mean over the window of 5 time points.

2 Model

The model consists of three states and a set of rules.

2.1 States

  • blue/brown (B)
  • pink (P)
  • red (R)

2.2 Rules

  • Simulation is carried out at a discrete time step of 1 min
  • There is a fixed time for nuclear envelope breakdown, \(t_0 = 0\); \(t_0\) can be changed if necessary
  • Time \(t_1\) is a random variable with exponential distribution with time scale \(\tau_1\)
  • Time \(\Delta t_2\) is a random variable with exponential distribution with time scale \(\tau_2\)
  • Time \(\Delta t_3\) is a random variable with exponential distribution with time scale \(\tau_3\)
  • The cell is in state B before time \(t_0-t_1\).
  • B\(\rightarrow\)P occurs after \(t_0-t_1\) with rate \(k_1\)
  • P\(\rightarrow\)R occurs after \(t_0-t_1 + \Delta t_2\) with rate \(k_2\)*
  • P\(\rightarrow\)B occurs after \(t_0-t_1 + \Delta t_3\) with rate \(k_3\)

* A switch, \(t_{2, ref}\) was introduced to the model, selecting the reference for time \(t_2\). In the default position (\(t_{2, ref} = 1\)) the B\(\rightarrow\)P activation time is \(t_0-t_1 + \Delta t_2\). When the switch is set to \(t_{2, ref} = 0\), the activation time is \(t_0 + \Delta t_2\), that is, it can only occur after the nuclear envelope breakdown. The idea of the switch was to separate pink and red curves, as observed in some data sets.

I use a Markov chain approach. The next state is generated from the current state based on rules outlined above. The rates, \(k\), are converted into probabilities over a given time step \(\Delta t\) as \(Pr = 1 - e^{1 - k\Delta t}\). The transition time, \(t_1\), is generated before the simulation starts for a given cell.

The cell timeline is repeated for \(n_{\rm cell}\) times and then the colour proportions are found at each time point.

2.3 Example

Here is an example of the model.

And here is an example of the the model where P\(\rightarrow\)R transition can occur only after \(t_0 = 0\).

2.4 Parameter tuner

There is a shiny app that allows tuning parameters in search for the best solution.

3 Fitting model to the data

I attempted fitting model to our data. It is not an easy task, as the model is stochastic. After some testing I found that modified BFGS (a quasi-Newton method, also known as a variable metric algorithm, by Broyden, Fletcher, Goldfarb and Shanno, 1970) which allows box constraints (Byrd et al. 1995), gives reasonable results. Again, due to stochasticity of the model, this algorithm often finds false local minima. I run it several times to find the best minimum. It is a crude and time-consuming method, but gives results better than manual tuning.

Fitting is constrained to time points between -50 and 30 min. The minimized quantity is

\({\rm rms} = \sqrt{\sum_{c \in \{B,P,R\}} \sum_i (O_{c,i} - E_{c,i})^2}\)

that is, the square root of the sum of squared residuals over all time points and all colours.

3.1 Default model

First, we fit the model with \(\tau_1\), \(\tau_2\), \(k_1\) and \(k_2\) free. To improve the chances of finding the correct global minimum, I run the fit with 10,000 cells and 100 tries.

3.1.1 Scramble

We suspect that \(t_0\) might not be exactly zero. Hence I fit the model with \(t_0 = 0, -5, -10, -15\) minutes.

I think the issue here is not \(t_0\) but the fact that the red curve growths is too fast and the model cannot do it.

3.1.2 All conditions

Here are fit results for all conditions with \(t_0\) fixed at zero and four parameters free: \(\tau_1\), \(\tau_2\), \(k_1\) and \(k_2\).

3.2 Switch: P\(\rightarrow\)R only after \(t_0\)

Here are results from a modified model, where P\(\rightarrow\)R transition can happen only after \(t_0\). By default \(t_0 = 0\) and the results are shown in the left panels below. As you can see, the red curve in the model lags behind the red curve in the data. Hence, I fitted the data again, but fixing \(t_0 = -10\) min this time (right panels). This aligns the red curves a little better.

We can also see that the pink curve is more peaked, as opposed to the smooth rise and decay in the default model.

4 Four-colour plots

Here I create cell-line plots for all data sets using four colours.

PDF files:

link
scramble
NCAPD2
NCAPD3
SMC2
WAPL24
WAPL48
MK1775
MK1775_ICRF193
RAD21
TT103
TT108

5 Brown density

Here I calculate the density of brown boxes and the distribution of intervals between them. If brown points are distributed randomly with density \(\lambda\), the intervals between them should follow the exponential distribution with PDF

\(f(d) = \lambda e^{-\lambda d}\)

where \(d\) is duration of the interval between consecutive brown points. Because some of the cells have much shorter data tracks, they will bias the interval distribution towards shorter intervals. The longer distances are missing from these short tracks. To account for this I made a simple model, in which I take cell tracks from the actual data, fill them randomly with brown points and calculate interval distribution. This is done 100 times to smooth the distribution.

The example below shows the result for scramble. The black bars represent data, the orange points show the predicted theoretical distribution \(f(d) = \lambda e^{-\lambda d}\) and the open blue bars show the simulated random distribution. As we can see, the simulated distribution predicts more short intervals with respect to the theoretical distribution, as expected.

I use a window of [-300, -20] min.

5.1 Brown-to-brown interval distributions

Below, are results for all conditions. For clarity, I don’t shows theoretical points. The number in the title is the brown point density \(\lambda\) (proportion of brown boxes).

5.2 K-S tests

Now, I do K-S test between each pair of distirbutions. P-values are in the table below.

scramble NCAPD2 NCAPD3 SMC2 WAPL24 WAPL48 MK1775 MK1775_ICRF193 RAD21 TT103 TT108
scramble
0.19 0.56 0.23 0.062 0.088 0.15 0.49 0.79 0.0017 0.12
NCAPD2 0.19
0.3 0.8 0.42 1 0.99 1 0.24 0.069 0.82
NCAPD3 0.56 0.3
0.24 0.27 0.5 0.39 0.31 0.66 0.017 0.52
SMC2 0.23 0.8 0.24
0.7 0.92 0.96 1 0.092 0.043 1
WAPL24 0.062 0.42 0.27 0.7
0.68 0.81 0.88 0.29 0.034 0.46
WAPL48 0.088 1 0.5 0.92 0.68
0.98 0.98 0.12 0.57 0.98
MK1775 0.15 0.99 0.39 0.96 0.81 0.98
1 0.13 0.76 1
MK1775_ICRF193 0.49 1 0.31 1 0.88 0.98 1
0.19 0.48 1
RAD21 0.79 0.24 0.66 0.092 0.29 0.12 0.13 0.19
0.00073 0.17
TT103 0.0017 0.069 0.017 0.043 0.034 0.57 0.76 0.48 0.00073
0.29
TT108 0.12 0.82 0.52 1 0.46 0.98 1 1 0.17 0.29

5.3 S and G2 data